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We study the origin of the cooperative nature of spin crossover (SC) between low spin (LS) and 
high spin (HS) states from the view point of elastic interactions among molecules. As the size of each 
molecule changes depending on its spin state, the elastic interaction among the lattice distortions 
provides the cooperative interaction of the spin states. We develop a simple model of SC with intra 
and intermolecular potentials which accounts for the elastic interaction including the effect of the 
inhomogeneity of the spin states, and apply constant temperature molecular dynamics based on the 
Nose-Hoover formalism. We demonstrate that, with increase of the strength of the intermolecular 
interactions, the temperature dependence of the HS component changes from a gradual crossover 
to a first-order transition. 

PACS numbers: 75.30.Wx 75.50.Xx 75.60.-d 64.60.-i 



The discovery of LIESST (light-induced excited spin 
state trapping) 1 -* 2 - phenomena has accelerated studies of 
functional spin-crossover (SC) molecular solids. SC com- 
pounds have been studied intensively not only because 
of their potential applicability to novel optical devices, 
e.g., optical data storage and optical sensors, etc, but 
also because of the fundamental scientific interest in the 
mechanism of the phase transition and the accompanied 
non-linear relaxation processes 1 -* 2 -^ 4 ^ 5 -^. 

To control electronic and magnetic states of SC com- 
pounds, it is important to understand the bistable nature 
of these compounds. The SC transition between the low- 
spin (LS) and high-spin (HS) states can be induced by 
change of temperature, pressure, magnetic field, light- 
irradiation, etc. It has been clarified that the interaction 
between spin states causes various types of cooperative 
phenomena between the LS and HS phases 8 . 

In order to take into account the cooperativity in the 
SC transition, Wajnflasz and Pick (WP) proposed an 
Ising model 9 , in which the spin state is described by a 
fictitious spin (<r = —1 (1) for the LS (HS) state) and 
the short-range interactions J between the spin states 
are introduced in the form TL = —J^( i j)ViVj J r — 
\k^T In g)<Ji, where g is the degeneracy ratio between the 
HS and LS states. Using this model, the change between 
a gradual crossover and a first-order transition has been 
well explained as a function of the parameters J, A and g. 
So far the WP model and its extensions called "Ising-like 
models" have been widely used for the description of the 
SC transition and related relaxation phenomena includ- 
ing photoinduced effects. Although the Ising-like models 
have captured several important feature a 10 ^ 1 ^ 2 ! 13 ! 14 , the 
origin of the parameters remains unclear due to the dras- 
tic simplifications involved. 

The importance of the elastic interaction in the SC 



transition has been investigate d 15 ! 16 iiLl&l&2&2L22 and 
the elastic energy of the system with the density distri- 
bution of the LS and HS sites has been phenomenolog- 
ically analyzed 19 . The dependences of elastic constants 
on the spin state have been also investigated in a one- 
dimensional (ID) two- level model 21 and in a ID vibronic 
coupling model 22 . However, local degrees of freedom 
(change of lattice) can be traced out in one dimension. 
Thus, no phase transition occurs in one dimension. 

Through the electron-distortion interaction, that is the 
vibronic coupling, the size of the molecule changes with 
the spin state. The distance (relative coordinate) be- 
tween the central transition metal and the surround- 
ing ligands changes. This distortion causes interactions 
among the spin states of molecules as depicted in Figs. Q] 
(a), (b), and (c). In the present study, we focus on the 
lattice distortions in higher dimensions caused by the dif- 
ference of the molecular sizes due to the different spin 
states. These local distortions interact with one another 
elastically which causes a long range effective interaction 
between the spin states. 

We perform molecular dynamics (MD) simulations on 
a 2D system with a simple square lattice. The in- 
tramolecular potential energy depending on the molec- 
ular size is given by a double-well adiabatic potential 
V^ mtra (ri), which is a function of the radius n of the zth 
molecule. Let pi be the corresponding relative momen- 
tum and let m be the reduced mass. 

We set an intermolecular binding interaction be- 
tween SC molecules (the zth and jth molecules) as 
ymter( x . x h n,rj), where X { = (Xi,Yj) is the co- 
ordinate of the center of the ith molecule (Fig. [2] (a)). 
The corresponding momentum is P{ = (Px^Pyi) and 
the mass of the molecule is M . 

To model this scenario we apply the following Hamil- 
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FIG. 1: (a) a schematic picture of the LS molecule and the HS 
molecule. The HS molecule is bigger than the LS molecule, 
(b) ((c)) shows a schematic picture of the lattice distortion 
for a HS (LS) molecule and surrounding LS (HS) molecules. 
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For simplicity we here consider only one symmetric 
vibration mode (isotropic volume expansion of molecules) 
as an active dominant mode2£. 

In order to clarify the effect of distortion, we adopt 
intermolecular binding potentials independent of the 
molecular states, although it is expected that the inter- 
molecular binding is looser for HS molecules than for LS 
molecules. 

As the intramolecular potential, we adopt a double well 
parabolic function V(x), where x is denned as the differ- 
ence of the radius from that of the ideal LS state. V(x) 
has minima at x = (ideal LS) and x = 1 (ideal HS). 
Setting tls = 9 and ths = 10 for the ideal radius of the 
LS molecule and that of the HS molecule, respectively, 
the radius of the molecule is r = tls + x. 

When a parabolic potential for the LS state (y = ax 2 ) 
and that for the HS state (y = b(x — c) 2 + d) are mixed 
by off-diagonal element J, the lowest potential function 
with coefficient A is given by 



V ( x ) = ~^{ d + K c ~ x ) 2 + ax<2 



(2) 



V4J 2 + (d + b(c - x) 2 - ax 2 ) 2 }. 



Because the entropy of the harmonic oscillator (H = 
£ + \Kx 2 ) is S = Nk B (l - ln^^/l), the entropy 
difference between the LS and HS states is given by 



AS = S HS - S LS = Nk B In 



-Khs 



(3) 
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FIG. 2: (a) a schematic picture of the model. (X, Y) is 
the coordinate of the center of each molecule and r is its 
radius, (b) Intramolecular potential, where x is the growth 
of r from tls- (c) Temperature dependence of (x) without 
intermolecular interactions. 



Thus the ratio of the degeneracy between the HS and 

LS states is g = yj = y^, which corresponds to g 

in the WP-model. We take A = 10, a = 10, b = 0.1, 
c = 1.0, d = 0.1, J = 0.04, which gives g = 10. In real- 
istic materials there are several sources of the difference 
of the entropy. However, in order to focus on only the 

lattice effect, we adopt a large value of y^J^f (= 10)- I n 

Fig 12(b), the intramolecular potential V(x) is depicted. 
The energy difference between the LS stable point and 
the HS stable point is A^ls-hs = 1.075 and the energy 
difference between the LS stable point and the unstable 
point (x c = 0.19) is AE acti = 1.654. The temperature 
dependence of the statistical average of x, i.e., 



(x) = 



Trxexp(-/?7i sys tem) 
Trexp(-/?W system ) 



(4) 



for V^ ter = 0, calculated by a numerical integration, is 
given in Fig. [2] (c). 

Next, we consider the intermolecular 
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(ith and jth molecules). We take 



y/f er ( Xi 



i r o)= f(dij), 



(5) 



where d^ = | Xi — Xj \ — (r^ + rj) . We treat phenomena 
in which the lattice distortion is not so large and it does 
not break the lattice structure. In this case, the quali- 
tative features of the phenomena do not depend on the 
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details of the potential form, and thus we adopt one of 
the simplest forms 



/(«)= D(e a >- 



-b'(u-uo) 



)■ 



(6) 



where a' = 0.5 and b' = 1.0. iio is a constant such that 
f(u) has the minimum at u = 0. When molecules (circles 
in Fig. [2] (a)) z and j contact each other (dij — 0), the 
function has a minimum value. 

In order to maintain the crystal structure (the coor- 
dination number), we introduce a potential between the 
next-nearest neighbors (i and /c, see Fig. [2] (a)) 



K fe nter ( Xi 



X k ,ri,r k )= f(d ik -Ar) 



(7) 



with a 1 = 0.1 and b' = 0.2, which is much smaller than 
that of the nearest neighbors. Next-nearest neighbors 
do not contact each other as depicted in Fig. [2] (a), and 
there is a spatial gap between them. For simplicity, we 
assume here that next-nearest neighbors are most stabi- 
lized when the gap is Ar = 2(v // 2 — l)f, where we take 
r = (^ls + ^hs)/2 although Ar can be temperature de- 
pendent. We focus on the dependence of the spin state on 
the strength of the intermolecular interaction, and thus 
we study the dependence on D. Common D is used for 
both Eqs. © and ®. 

To study the temperature dependence, we adopt the 
Nose-Hoover method 24,25 to generate the canonical en- 
semble for a given temperature T. The Hamiltonian of 
the thermal reservoir is given by 
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(8) 



where s is a scaling factor, P s is the conjugated momen- 
tum of s and Q is an effective mass associated with s. 
Therefore, the total Hamiltonian including the effect of 
thermal reservoir is given by H tota i = H syst em + Htherm- 
Applying the Nose-Hoover formalism to the present 
system, the time evolution of the system is realized ac- 
cording to the following equations of motion. 
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where y mter stands for the summation of the intermolec- 
ular potentials for the nearest and next-nearest pairs, and 
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FIG. 3: Temperature dependences of (x) for the values of (a) 
D = 10, (b) D = 20, (c) D = 28, and (d) D = 42. The open 
red circles (blue squares) denote (x) in the warming (cooling) 
process. 



We adopt x as a parameter to characterize the spin 
state. We study the open-boundary system of L 2 — 26 x 
26 molecules. We warm up the system from T = 0.1 to 
2.0 in steps of increment 0.1, and cool it down to T = 
0.1. At each temperature, 40000 MD steps are discarded 
as transient time and subsequent 20000 MD steps are 
used to measure x with the time step At = 0.01. We 
employ an operator decomposition method in which the 
numerical error is of the order 0(At 3 ). We set m — 1.0, 
M = 1.0, and Q = 1.0. ((x) does not depend on m, M, 
and Q in the equilibrium state.) 

In Figs. [3] (a), (b), (c) and (d), the temperature depen- 
dences of (x) are shown. When D = 10 (Fig. [3] (a)), (x) 
in the warming process and that in the cooling process 
overlap, indicating that a smooth (gradual) SC crossover 
is realized. When the interaction parameter becomes 
larger: D = 20 (Fig. [3] (b)), variation of (x) becomes 
sharper, which implies that the SC transition becomes 
more cooperative. 

When D = 28 (Fig. [3] (c)), a clear hysteresis loop of 
(x) is found. As the interaction parameter increases fur- 
ther: D = 42 (Fig. [3] (d)), the hysteresis width becomes 
larger. Here, we found that when the interaction between 
molecules becomes large, the SC transition changes from 
a gradual crossover to a first order transition. The criti- 
cal Value of D is Identical — 20. 

In Figs. [H (a) and (b), snapshots of the complete LS 
state and the complete HS state are shown, in which the 
system length changes by 11 %. In Figs. [H (c) and (d), 
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FIG. 4: Snapshots of configurations. HS molecules (red cir- 
cles) are allocated when r is larger than r = tls + x c . LS 
molecules are drawn by blue circles, (a) Complete LS phase, 
(b) Complete HS phase, (c) A snapshot of configuration at 
T = 0.6 in the system of D = 10 (L 2 = 16 2 ), where 80 
molecules are in the HS state, (d) A snapshot of configura- 
tion at T = 1.0 in the system of D = 42 (L 2 = 16 2 ), where 
79 molecules are in the HS state. 

a snapshot of configuration at T — 0.6 for the parameter 
D = 10 and that at T = 1.0 for D = 42 are given, where 
the concentration of HS molecules is about 30% in both 



configurations. Although the number of HS molecules is 
almost the same, the average cluster size of HS in (d) 
is bigger than in (c). This indicates that there is higher 
correlation between spin states of molecules in the case of 
strong intermolecular interaction (case (d)), which pro- 
motes first-order transition. 

In this study, we investigated the cooperativity of spin- 
crossover phenomena induced by the elastic interaction 
among lattice distortions which are triggered by the dif- 
ference of molecular sizes caused by the different spin 
states. This effect is inherent to the high dimensional- 
ity (2D and 3D). The present 2D model can be applied 
straightforwardly to the 3D case. Although the lattice 
relaxation through a change of molecular sizes has been 
studied phenomenologically by a mean- field treatment 19 , 
as far as we know, this is the first attempt to investigate 
the cooperativity attributed to the effect of local distor- 
tions (fluctuation) and that of the propagation to the 
overall lattice. 

The present work was supported by Grant-in- Aid for 
Scientific Research on Priority Areas and for Young Sci- 
entists (B) from MEXT of Japan, and supported by 
Grant-in-Aid from Ministere de l'Education Nationale 
and CNRS (PICS France- Japan Program) of France. The 
present work was also supported by the MST Founda- 
tion. The numerical calculations were supported by the 
supercomputer center of ISSP of Tokyo University. 



* Corresponding author. Email address: 

nishino . masamichi@ nims . go . j p 

1 P. Giitlich, A. Hauser and H. Spiering, Angew. Chem. Int. 
Ed. Engl. 33, 2024 (1994) and references therein. 

2 S. Decurtins, P. Giitlich, C.P.Kohler, H. Spiering, and A. 
Hauser, Chem. Phys. Lett. 105, 1 (1984). 

3 O. Kahn and C. Jay Martinez, Science 279, 44 (1998). 

4 J. F. Letard, P. Guionneau, L. Rabardel, J. A. K. Howard, 
A. E. Goeta, D. Chasseau, and O. Kahn, Inorg. Chem., 
37, 4432 (1998). 

5 A. Hauser, J. Jeftic, H. Romstedt, R. Hinek and H. Spier- 
ing, Coord. Chem. Rev. 190-192, 471 (1999). 

6 J. A. Real, H. Bolvin, A. Bousseksou, A. Dworkin, O. 
Kahn, F. Varret, and J. Zarembowitch, J. Am. Chem. Soc. 
114, 4650 (1992). 

7 N. Shimamoto, S. Ohkoshi, O. Sato, and K. Hashimoto, 
Inorg. Chem. 41, 678 (2002). 

8 M. Sorai and S. Seki, J. Phys. Chem. Solids, 35, 555 (1974). 

9 J. Wajnflasz and R. Pick, J. Phys. (Paris). Colloq 32, Cl- 
91 (1971). 

10 A. Bousseksou, J. Nasser, J. Linares, K. Boukheddaden, 
and F. Varret, J. Phys. I France 2, 1381 (1992); A. 
Bousseksou, F. Varret, and J. Nasser, ibid. 3, 1463 (1993); 
A. Bousseksou, J. Constant-Machado, and F. Varret, ibid. 
5, 747 (1995). 

11 K. Boukheddaden, I. Shteto, B. Hoo, and F. Varret., Phys. 
Rev. B 62, 14796 (2000); ibid. 14806 (2000). 



12 M. Nishino, S. Miyashita, and K. Boukheddaden, J. Chem. 
Phys. 118, 4594 (2003). 

13 M. Nishino, K. Boukheddaden, S. Miyashita, and 
F.Varret., Phys. Rev. B 72, 064452 (2005). 

14 S. Miyashita, Y. Konishi, H. Tokoro, M. Nishino, K. 
Boukheddaden, and F.Varret, Prog. Theor. Phys. 114, 719 
(2005). 

15 R. Zimmermann and E. Konig, J. Phys. Chem. Solids 38, 
779 (1977). 

16 T. Kambara, J. Chem. Phys. 70, 4199 (1979). 

17 P. Adler, L. Wiehl, E. MeiBner, CP Kohler, H. Spiering, 
and P. Giitlich, J. Phys. Chem. Solids 48 517 (1987). 

18 N. Willenbacher and H. Spiering. J. Phys. C: Solid State 
Phys. 21, 1423 (1988). 

19 A. L. Tchougreeff and M. B. Darkhovskii, Int. J. Quant. 
Chem. 57, 903 (1996). 

20 H. Spiering, K. Boukheddaden, J. Linares, and F. Varret, 
Phys. Rev. B. 70, 184106 (2004). 

21 J. A. Nasser, K. Boukheddaden, and J. Linares. Eur. Phys. 
J. B. 39, 219 (2004). 

22 K. Boukheddaden, Prog. Theor. Phys. 112, 205 (2004). 

23 A. Bousseksou, J. J. McGarvey, F. Varret, J.A. Real, J.-P. 
Tuchagues, A.C. Denis, and M. L. Boillot, Chem. Phys. 
Lett. 318, 409 (2000). 

24 S. Nose, J. Chem. Phys. 81, 511 (1984). 

25 W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 



